RANK signaling in osteoclast precursors results in a more permissive epigenetic landscape and sexually divergent patterns of gene expression

Background Sex is an important risk factor in the development of osteoporosis and other bone loss disorders, with women often demonstrating greater susceptibility than men. While variation in sex steroids, such as estradiol, accounts for much of the risk, there are likely additional non-endocrine factors at transcriptional and epigenetic levels that result in a higher rate of bone loss in women. Identification of these factors could improve risk assessment and therapies to preserve and improve bone health. Methods Osteoclast precursors were isolated male and female C57Bl/6 mice and cultured with either MCSF alone or MCSF and RANKL. Following the culture period RNA was isolated for RNA sequencing and DNA was isolated for tagmentation and ATAC sequencing. RNA-Seq and ATAC-seq were evaluated via pathway analysis to identify sex- and RANKL-differential transcription and chromatin accessibility. Results Osteoclasts demonstrated significant alterations in gene expression compared to macrophages with both shared and differential pathways between the sexes. Transcriptional pathways differentially regulated between male and female cells were associated with immunological functions with evidence of greater sensitivity in male macrophages and female osteoclasts. ATAC-Seq revealed a large increase in chromatin accessibility following RANKL treatment with few alterations attributable to sex. Comparison of RNA-Seq and ATAC-seq data revealed few common pathways suggesting that many of the transcriptional changes of osteoclastogenesis occur independently of chromatin remodeling.


INTRODUCTION
Osteoclasts are multinucleated hematopoietic-lineage cells responsible for the dissolution of bone matrix required for both physiological bone growth and remodeling (Teitelbaum, 2007). While indispensable for normal bone homeostasis, excessive osteoclastic resorption can lead to insufficient bone mineral density, joint destruction, tooth loss, and fracture (Ensrud & Crandall, 2017). Multiple genetic and environmental factors contribute to overall risk of pathologically elevated osteoclast activity, but the greatest contributors are advanced age and female sex (Watts et al., 2008;Johannes de Villiers, 2009). Osteoporosis, a generalized decrease in bone mass to at least 2.5 standard deviations below that of an average young adult, results in a lifetime risk of fracture of one in five among men over fifty and one in three among women over fifty (De castro machado, Hannon & Eastell, 1999;Yelin, Weinstein & King, 2016;Sebbag et al., 2019). The high rate of osteoporosis resulting in both direct medical costs and loss of productivity and quality of life is an impetus for improving understanding of additional risk factors and identifying new therapeutic targets.
The disproportionate osteoporosis risk among women is explained mostly by the regulatory role of estrogens in osteoclast differentiation and function. Estradiol limits resorption by both reducing differentiation and limiting the lifespan of osteoclasts (Shevde et al., 2000;Krum et al., 2008). Prior to menopause, ovary-derived estradiol protects young female bone mass. After menopause, which is marked by cessation of ovarian follicle development and a precipitous decline in circulating estradiol, osteoclasts can become disinhibited. If the hyperactivity of postmenopausal osteoclasts cannot be overcome by increased bone-forming activity of osteoblasts, progressive bone loss and osteoporosis follows (Recker et al., 2004;Iki et al., 2004).
Men experience osteoporosis risk differently. Although circulating estradiol level in men of all ages is comparable to that of post-menopausal women, aromatase enzyme (encoded by the CYP19A1 gene) present in the bone converts testosterone to estradiol locally, where it exerts its protective effect (Sasano et al., 1997;Gennari, Nuti & Bilezikian, 2004;Khosla et al., 2005aKhosla et al., , 2005b. In contrast to women and estradiol, men do not experience a rapid decrease in testosterone at an age-dependent milestone like menopause, but rather demonstrate increasing variation beginning in the second decade of life (Kelsey et al., 2014). Sex-specific endocrine differences likely account for much of sex-dependent osteoporosis risk, with bone parameters of men in their eighth decade approaching parity with women in their sixth (Alswat, 2017).
While estrogen deficiency explains much of the sex-differential osteoporosis risk, the regulators of osteoclast differentiation and activity are diverse, and there are likely sexually divergent contributors independent of sex steroids. Osteoclasts differentiate from macrophage-like precursors via stimulation with Receptor Activator of Nuclear Factor κB Ligand (RANKL). Upon binding of RANKL to its cognate receptor, RANK, Tumor necrosis factor Receptor Associated Factor (TRAF)-dependent signaling cascades converge on multiple transcription factors including JNK, ERK, p38, and NFκB (Feng, 2005). This receptor-proximal signaling alters expression of multiple genes; chief among them is upregulation of Nuclear Factor of Activated T Cells 1 (NFATc1), which drives much of the osteoclast differentiation program (Kim & Kim, 2014, p. 1). Expression of NFATc1 and other osteoclastic genes is also subject to epigenetic regulation (Rohatgi et al., 2018;Das et al., 2018;Shin et al., 2019;Astleford et al., 2020). The complexity of osteoclastic gene expression presents opportunities for variation in degree of differentiation and level of resultant resorptive activity, and multiple pathways other than RANKL/RANK can modulate osteoclastogenesis (Liu et al., 2009;Goel et al., 2019). From this we hypothesized that there are sexually divergent patterns of gene regulation in the differentiation of osteoclasts.
Here we report the findings of our analyses of gene expression and chromatin accessibility in bone marrow macrophages and osteoclasts extracted from female and male mice.

MATERIALS AND METHODS
Osteoclast precursor culture Mice used as cell sources were euthanized and provided by the staff of the Eastern Washington University vivarium (PHS animal welfare assurance number D19-01059). Primary mouse bone marrow cells were obtained by flushing a-MEM medium (Minimum Essential Medium Eagle (M0894;MilliporeSigma,St. Louis,MO,USA), 20 mM GlutaMAX (35050061; Thermo Fisher, Waltham, MA, USA), and 10% heat-inactivated fetal bovine serum (F-500-D; Atlas Biologicals, Fort Collins, CO, USA)) through the femurs and tibias of 3-month-old C57Bl6 female and male mice (three each for RNA-sequencing; two each for ATAC-Sequencing) (Kaur et al., 2017). Bone marrow cells were cultured at 37 C and 5% CO2 overnight in tissue culture-treated dishes to allow adherent cells to attach. The next day, non-adherent cells were transferred to non-treated suspension culture dishes and supplemented with 25 ng/mL Macrophage Colony-Stimulating Factor/MCSF (200-08; Shenandoah Biotechnology, Warminster, PA, USA) to induce the differentiation of macrophages, which can attach even to non-treated culture surfaces. After 24 h, the MCSF-containing medium was refreshed to remove non-macrophage cells, and macrophages were maintained for an additional 48 h to allow them to proliferate prior to experimental treatments.

RNA-Sequencing (RNA-Seq)
Appropriate sample sizes were calculated using RNASeqPower with the following parameters: average read depth = 20 million, biological coefficient of variation = 0.1, effect size = 2, alpha = 0.05, power = 0.9 (Hart et al., 2013;Therneau & Stephen, 2022). From this, a sample size of three measurements per group were determined as sufficient. Osteoclast precursors from male and female mice were divided into two treatment groups: one group of cells (macrophages) were maintained in a -MEM with 25 ng/mL MCSF for 72 h; the other group of cells (osteoclasts), were maintained in a -MEM with 25 ng/mL MCSF (200-08; Shenandoah Biotechnology, Warminster, PA, USA) and 100 ng/mL RANKL (200-04; Shenandoah Biotechnology, Warminster, PA, USA) for 72 h. Media were refreshed at 48 h. At the conclusion of the culture period, cells were lysed in TRI reagent and RNA was extracted using the Direct-zol RNA miniprep kit (R2051; Zymo Research, Irvine, CA, USA). RNA was extracted from each group (female-macrophages, female-osteoclasts, male-macrophages, and male-osteoclasts) in triplicate. RNA samples were transferred to GeneWiz/Azenta for cDNA synthesis, indexing, sequencing, and basic RNA-seq analysis (de-multiplexing, alignment, transcript identification, and differential expression analysis). Raw sequence reads were trimmed to remove adaptors and poor sequence quality nucleotides using Trimmomatic v.0.36, trimmed reads were mapped to the GRCm38 Mus musculus reference genome using STAR aligner v.2.5.2b, unique gene hit counts from reads falling within exon regions were calculated with the featureCounts module of Subread package v.1.5.2, and differential expressions were calculated using DESeq2 with pvalues and log2 fold changes determined via the Wald test. Differential expression data is included as Supplemental Files. RNA sequencing data is available at NCBI GEO, accession number GSE216929

Assay for transposase-accessible chromatin-sequencing
Osteoclast precursors from male and female mice were divided into two groups: one group of cells (MCSF only) were maintained in a-MEM with 25 ng/mL MCSF for 24 h; the other group of cells (MCSF+RANKL),were maintained in a-MEM with 25 ng/mL MCSF and 100 ng/mL RANKL for 24 h. At the conclusion of the culture period, nuclei were isolated and used in tagmentation and indexing reactions using the ATAC-Seq Kit (53150, Active Motif, Carlsbad, CA, USA). ATAC-seq libraries were prepared from each sex and treatment group (female-MCSF, female-MCSF+RANKL, male-MCSF, male-MCSF +RANKL) in duplicate. ATAC-seq libraries were transferred to GeneWiz/Azenta, which performed the sequencing reactions and basic ATAC-seq analysis (de-multiplexing, alignment, peak calling, and peak differential analysis). Sequencing adaptors and low-quality nucleotides were removed with Trimmomatic v.0.38, and reads were aligned to the Mus musculus mm10 reference genome using bowtie2. Aligned reads were filtered with samtools v1.9 to preserve alignments with a minimum mapping quality of 30, are aligned concordantly, and are the primary called alignments; PCR or digital duplicates were marked with Picard v2.18.26 and removed. Mitochondrial reads and reads mapped to unplaced contigs were also removed. Open chromatin regions were identified through peak calling with MACS2 v2.1.2, and consensus peaks from sample duplicates were merged and kept for downstream analysis. After reads falling beneath peaks were counted, differential peak calling (corresponding to differentially accessible chromatin regions) was performed using the R package Diffbind. Differential peak call data are provided as Supplemental Files.

Pathway analyses of differential sequencing data
Differential data with non-adjusted p-values were uploaded to iPathwayGuide (Advaita Bioinformatics; http://www.advaitabio.com/ipathwayguide). iPathwayGuide is a web-based bioinformatics platform that applies Impact Analysis to differential datasets, which considers both over-representation of differential genes in a particular pathway and computed perturbation to generate pathway specific p-values (Tarca et al., 2009;Ahsan & Drăghici, 2017). p-Values for individual genes were adjusted via False Discovery Rate. Venn diagrams were generated with eulerr.co (Larsson, Godfrey & Gustafsson, 2021).  Similarly, when comparing cells of the same differentiation state, sex of the cells was the principal component ( Fig. 2B).

Gene expression patterns in macrophages and osteoclasts cluster by sex and differentiation state
RNA-sequencing of macrophages and osteoclasts identifies expected differential gene expression Figure 3 depicts biclustering analyses of 30 differentially expressed genes with the lowest adjusted p-values according to differentiation state and sex. Among both male and female cells, matrix metalloproteinase 9 (Mmp9), cathepsin K (Ctsk), and tartrate-resistant acid phosphatase (Acp5), which are involved in osteoclast function, were among the most significantly differential (Fig. 3A). When comparing cells of the same differentiation state and opposite sexes, Ubiquitously Transcribed Tetratricopeptide Repeat Containing, Y-Linked (Uty), DEAD-Box Helicase 3 Y-Linked (Ddx3y), and Eukaryotic translation initiation factor 2 subunit 3, Y-linked (Eif2s3y), which are restricted to the Y chromosome, were found to be more highly expressed in male cells. Conversely, X-inactive specific transcript (Xist), which is expressed only in cells with more than one X chromosome, was more highly expressed in female cells (Fig. 3B). The identification of expected differences in gene expression by the unbiased RNA-seq method increases confidence in the method and subsequent analyses to identify novel differential gene expression based on sex and differentiation state.
Sex-independent and sex-dependent differential pathways in osteoclastogenesis Using a p-value threshold of 0.001 and a minimum fold change of 1.5, 2,280 differentially expressed genes were identified in female cells and 2,844 in male (Fig. 4A). Following false discovery rate correction, pathway analysis of differential gene expression according to differentiation state identified 22 altered pathways common to both female and male cells, 45 specific to male cells, and two specific to female cells (Fig. 4B). These pathways are listed in Table 1.
Differentiation-independent and differentiation-dependent differential pathways in female and male cells Using a p-value threshold of 0.05 and a minimum fold change of 0.6, 1,347 differentially expressed genes were identified in macrophages and 1,618 in osteoclasts (Fig. 5A). Following false discovery rate correction, pathway analysis of differential gene expression according to sex identified 32 altered pathways common to both macrophages and osteoclasts, 15 specific to macrophages, and 34 specific to osteoclasts (Fig. 5B). These pathways are listed in Table 2.
Differential transposase accessible chromatin regions regulated by RANKL in female and male cells ATAC-seq revealed that most differences in transposase-accessible chromatin were in promotor and distal intergenic regions (Fig. 6). Using a p-value threshold of 0.01 and a minimum accessibility fold difference of 2.7 (female) or 2.4 (male), there were 4,957 genes with differentially accessible chromatin in female MCSF+RANKL-treated cells and 4,994 in male (Fig. 7A). Most differentially accessible genes were more accessible following MCSF+RANKL treatment with more accessible genes accounting for 95% in female cells and 94% in male cells. Without FDR correction, there were 17 differentially accessible pathways specific to female cells, six specific to male cells, and five common to both sexes (Fig. 7B). These pathways are listed in Table 3.

Differential gene expression and gene accessibility during osteoclastogenesis
Meta-analysis of macrophage vs osteoclast RNA-seq and MCSF vs MCSF+RANKL ATAC-seq data revealed 51 gene expression specific pathways, 17 gene accessibility specific pathways, and five pathways common to both in female cells (Fig. 9A). These pathways are listed in Table 5. Among male cells, there were 89 gene expression specific pathways, seven gene accessibility specific pathways, and four common pathways (Fig. 9B). These pathways are listed in Table 6.

DISCUSSION
In this study, we applied high-throughput sequencing and bioinformatic analyses to characterize the transcriptional profile and epigenetic landscape of osteoclastogenesis in cells derived from female and male mice. Throughout this study, we focused our analyses on differences related to sex and differentiation state with a goal of identifying novel genes and pathways that may influence osteoclast differentiation in a sexually divergent manner. Our initial analyses correctly identified well-established genes involved in osteoclast function as well as known sexually divergent genes such as XIST (overexpressed in XX cells) and Uty (found only on Y chromosomes). That our unbiased analyses correctly identified anticipated variation at both the transcriptional and epigenetic level supports their ability to identify novel differences as well.
To manage the large datasets produced in this study, we utilized iPathwayGuide, which applies Impact Analysis to assign genes to pathways described by the Kyoto Encyclopedia of Genes and Genomes (KEGG) and statistically evaluates pathways according to fold changes and p-values of differentially expressed genes (Tarca et al., 2009;Ahsan & Drăghici, 2017). We deliberately selected low fold change thresholds for differentially expressed genes when performing pathway analyses for two reasons: first, there is evidence that small fold changes in gene expression can result in significant biological impacts and setting fold change thresholds too high can lead to disposal of potentially relevant genes, and second, by providing more genes to iPathwayGuide, the analysis was better equipped  to evaluate pathway-level changes (St. Laurent et al., 2013). The top two sex-independent pathways altered by osteoclast differentiation, ranked by FDR-corrected p-value, were metabolic pathways and axon guidance. With respect to metabolic pathways, alterations were diverse with most changes falling within pathways of fatty acid metabolism, citric acid  cycle, and oxidative phosphorylation. Almost without exception, when gene expression was altered in these pathways, it was increased. This concurs with prior findings of enhanced ATP production and mitochondrial activity in mature osteoclasts (Li et al., 2020;Da, Tao & Zhu, 2021). Some members of the axon guidance pathway have already been implicated in bone homeostasis such as the Eph/Ephrin and Slit/Robo axes (Matsuo & Otaki, 2012;Kim et al., 2018, p. 3;Jiang, Sun & Huang, 2022). In addition to identifying increased EphA, EphrinA, and EphrinB gene expression and decreased Slit1 expression in  Table 3. Full-size  DOI: 10.7717/peerj.14814/ fig-7 osteoclasts, we found altered expression of other axon guidance membrane receptors including Frizzled3 (Fzd3) and biregional cell adhesion molecule-related/down-regulated by oncogenes binding protein (Boc), which have not been directly interrogated for roles in osteoclast function. Alterations in pathways associated with immune function was a common finding when comparing male and female-derived cells, which is in agreement with prior published findings (Gal-Oz et al., 2019;Mun et al., 2021). Interestingly, for multiple pathways linked to immune responses to diverse pathogens, expression of pattern recognition receptors is higher in male macrophages, but, after differentiation, female osteoclasts demonstrate higher expression. In the Toll-like receptor (TLR) pathway, for example, expression of all TLRs except TLR5 (for which there was no significant difference) was higher in female osteoclasts. By contrast, expression of TLR1, TLR2, and TLR9 was higher in male macrophages. These findings suggest that male naïve macrophages are more sensitive to TLR ligands and agree with prior findings of sexually dimorphic responses to inflammatory signals and induction of TLR expression in conjunction with pro-osteoclastogenic pathways (Marriott, Bost & Huet-Hudson, 2006;Li et al., 2009;Barcena et al., 2021;Chen, Lainez & Coss, 2021;Song et al., 2022). The potency of some TLR ligands, such as bacterial lipopolysaccharide (LPS), to stimulate osteoclastogenesis of osteoclast precursors pre-committed with RANKL is well established (Liu et al., 2009). It has also been shown that there is a sexually divergent response to LPS by differentiating osteoclasts with female cells demonstrating more robust differentiation (Mun et al., 2021).
In addition to pathways with previously established roles in osteoclast function, our analyses identified additional pathways that have not been studied directly in osteoclasts.  Table 4.
Full-size  DOI: 10.7717/peerj.14814/ fig-8 For example, the Synaptic vesicle cycle (KEGG: 04721) was identified as a sex-independent up-regulated pathway in osteoclastogenesis. While this pathway has also appeared in analyses of differentially expressed genes in post-menopausal osteoporosis, roles of specific differential genes, such as Cplx2, Cacna1b, and Dnm1, have not been investigated (Zhu et al., 2018). This pathway and others not currently associated with macrophage or osteoclast function might represent the richest areas for discovery of novel gene functions, and we intend to investigate candidate genes of these pathways for deeper insights into osteoclast regulation.
As did the RNA-seq analyses, the ATAC-Seq studies identified expected sexually divergent chromosome accessibility with Y chromosome-specific loci overrepresented in male samples and X chromosome-specific loci overrepresented in female samples. The degree to which RANK signaling loosens chromatin is remarkable, with more accessible genes accounting for greater than 94% of the differences in both sexes. Given the balance of genes that are up-and down-regulated by RANK signaling, a similar pattern of chromatin accessibility might be expected. On the contrary, RANK signaling results in a more transcriptionally accessible epigenetic landscape. This might explain the ability of diverse inflammatory signals to complete osteoclastogenesis in RANKL-committed precursors-RANK signaling might render previously inaccessible pro-osteoclast genes available to transcription factors downstream from inflammatory receptors. By comparison, sex appears to have little influence in chromatin accessibility, suggesting that alterations in gene expression between male and female osteoclast precursors are driven primarily by transcriptional, rather than epigenetic, mechanisms.
Transcription appears to be a greater driver of osteoclastogenesis than epigenetic changes. In both female and male osteoclast precursors, pathway analysis revealed a greater number of significantly altered pathways following RNA-seq than ATAC-seq with minimal overlap between the two. Axon guidance was the only pathway with concordance between RNA-seq and ATAC-seq in both sexes, and while RNA-seq revealed a mixture of up-and down-regulated genes, ATAC-seq demonstrated an almost uniform increase in accessibility. While increased transcription in RANKL-treated cells may depend on increased chromatin accessibility, our data suggest that decreased transcription is not reliant on decreased chromatin accessibility.  Table 5. (B) RNA-Seq-specific, ATAC-specific, and common pathways in male cells. Pathways are listed in Table 6.

CONCLUSIONS
In this study we hypothesized that gene expression patterns and chromatin accessibility in osteoclast precursors are regulated by both RANK signaling and sex-specific factors. Using RNA-seq, ATAC-seq, and pathway analyses, we found that RANKL produces different transcriptional patterns in male-and female-derived cells, and these patterns suggest that male and female osteoclast precursors may respond differently to inflammatory signals. Furthermore, we found that RANK signaling results in a more permissive epigenetic landscape, with most loci becoming more accessible following treatment with RANKL.
Multiple studies have demonstrated that inflammatory signals including LPS, poly I:C, and TNF can support the differentiation of committed osteoclast precursors while those same factors prevent osteoclastogenesis of naïve precursors. Our data suggest that male and female osteoclast precursors may have different sensitivities to these and other factors. Future studies should investigate these potential differences, which may further explain sexually divergent risk of bone loss.